lambda = -1e6;
eta = 1.5; 
k = 0.1; % Step size
u = @(t) exp(lambda*t)*(eta - 1) + cos(t); % Exact solution
f = @(u,t) lambda*(u - cos(t)) - sin(t); % ODE definition

% Exact solution for plotting
x = linspace(0, 3, 1000);
y = u(x);

% Backward Euler
t_be = 0:k:3;
u_be = zeros(size(t_be));
u_be(1) = eta; % Initial condition

for n = 1:length(t_be)-1
    % Define the implicit equation to solve: F(u_{n+1}) = 0
    F = @(u_star) u_star - u_be(n) - k * f(u_star, t_be(n+1));
    
    % Solve using MATLAB's fsolve (requires Optimization Toolbox)
    options = optimset('Display','off');
    u_be(n+1) = fsolve(F, u_be(n), options); % Initial guess: u_be(n)
end

% Second-order trapezoidal
t_trap = 0:k:3;
u_trap = zeros(size(t_trap));
u_trap(1) = eta; % Initial condition

for n = 1:length(t_trap)-1
    % Define the implicit equation to solve: F(u_{n+1}) = 0
    F = @(u_star) u_star - u_trap(n) - k/2 * (f(u_star, t_trap(n+1)) + f(u_trap(n), t_trap(n)));
    
    % Solve using MATLAB's fsolve (requires Optimization Toolbox)
    options = optimset('Display','off');
    u_trap(n+1) = fsolve(F, u_trap(n), options); % Initial guess: u_be(n)
end


% Plot
figure; hold on;
plot(x, y, 'k', 'LineWidth', 1.5); % True solution
plot(t_be, u_be, 'b--', 'LineWidth', 1.5); % Backward Euler
plot(t_be, u_be, 'bs--', 'LineWidth', 1, 'MarkerSize', 5, 'MarkerFaceColor', 'b'); % 后向欧拉（蓝色方块虚线）
legend('True Solution', 'Backward Euler', 'Location', 'best');
title('Stiff ODE: Backward Euler vs Exact Solution');
xlabel('t'); ylabel('u');
grid on;


% Plot
figure; hold on;
plot(x, y, 'k', 'LineWidth', 1.5); % True solution
plot(t_trap, u_trap, 'r--', 'LineWidth', 1.5); % Second-order trapezoidal
plot(t_trap, u_trap, 'rs--', 'LineWidth', 1, 'MarkerSize', 5, 'MarkerFaceColor', 'r'); % Second-order trapezoidal
legend('True Solution', 'Second-order trapezoidal', 'Location', 'best');
title('Stiff ODE: Second-order trapezoidal vs Exact Solution');
xlabel('t'); ylabel('u');
grid on;